home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / newmat.exe / NEWMAT7.CXX < prev    next >
C/C++ Source or Header  |  1991-07-30  |  11KB  |  348 lines
  1. //$$ newmat7.cxx     Invert, solve, binary operations
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #include "include.hxx"
  6.  
  7. #include "boolean.hxx"
  8.  
  9. typedef double real;                 // all floating point double
  10.  
  11. #include "newmat.hxx"
  12.  
  13. #include "newmatrc.hxx"
  14.  
  15. //#define REPORT { static ExeCounter ExeCount(__LINE__,7); ExeCount++; }
  16.  
  17. #define REPORT {}
  18.  
  19.  
  20. /***************************** solve routines ******************************/
  21.  
  22. GeneralMatrix* GeneralMatrix::MakeSolver()
  23. {
  24.    REPORT
  25.    GeneralMatrix* gm = new CroutMatrix(*this);
  26.    MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
  27. }
  28.  
  29. GeneralMatrix* Matrix::MakeSolver()
  30. {
  31.    REPORT
  32.    GeneralMatrix* gm = new CroutMatrix(*this);
  33.    MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
  34. }
  35.  
  36. void CroutMatrix::Solver(MatrixRowCol& mcout, const MatrixRowCol& mcin)
  37. {
  38.    REPORT
  39.    real* el = mcin.store; int i = mcin.skip;
  40.    while (i--) *el++ = 0.0;
  41.    el += mcin.storage; i = nrows - mcin.skip - mcin.storage;
  42.    while (i--) *el++ = 0.0;
  43.    lubksb(mcin.store, mcout.skip);
  44. }
  45.  
  46.  
  47. // Do we need check for entirely zero output?
  48.  
  49. void UpperTriangularMatrix::Solver(MatrixRowCol& mcout,
  50.    const MatrixRowCol& mcin)
  51. {
  52.    REPORT
  53.    real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
  54.    while (i-- > 0) *elx++ = 0.0;
  55.    int nr = mcin.skip+mcin.storage; elx = mcin.store+nr; real* el = elx;
  56.    int j = mcout.skip+mcout.storage-nr; int nc = ncols-nr; i = nr-mcout.skip;
  57.    while (j-- > 0) *elx++ = 0.0;
  58.    real* Ael = store + (nr*(2*ncols-nr+1))/2; j = 0;
  59.    while (i-- > 0)
  60.    {
  61.       elx = el; real sum = 0.0; int jx = j++; Ael -= nc;
  62.       while (jx--) sum += *(--Ael) * *(--elx);
  63.       elx--; *elx = (*elx - sum) / *(--Ael);
  64.    }
  65. }
  66.  
  67. void LowerTriangularMatrix::Solver(MatrixRowCol& mcout,
  68.    const MatrixRowCol& mcin)
  69. {
  70.    REPORT
  71.    real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
  72.    while (i-- > 0) *elx++ = 0.0;
  73.    int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.store+i;
  74.    int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
  75.    while (j-- > 0) *elx++ = 0.0;
  76.    real* el = mcin.store+nc; real* Ael = store + (nc*(nc+1))/2; j = 0;
  77.    while (i-- > 0)
  78.    {
  79.       elx = el; real sum = 0.0; int jx = j++; Ael += nc;
  80.       while (jx--) sum += *Ael++ * *elx++;
  81.       *elx = (*elx - sum) / *Ael++;
  82.    }
  83. }
  84.  
  85. /******************* carry out binary operations *************************/
  86.  
  87. static void Add(GeneralMatrix*,GeneralMatrix*, GeneralMatrix*);
  88. static void Add(GeneralMatrix*,GeneralMatrix*);
  89.  
  90. static void Subtract(GeneralMatrix*,GeneralMatrix*, GeneralMatrix*);
  91. static void Subtract(GeneralMatrix*,GeneralMatrix*);
  92. static void ReverseSubtract(GeneralMatrix*,GeneralMatrix*);
  93.  
  94. static GeneralMatrix* GeneralAdd(GeneralMatrix*,GeneralMatrix*,MatrixType);
  95. static GeneralMatrix* GeneralSub(GeneralMatrix*,GeneralMatrix*,MatrixType);
  96. static GeneralMatrix* GeneralMult(GeneralMatrix*,GeneralMatrix*,MatrixType);
  97. static GeneralMatrix* GeneralSolv(GeneralMatrix*,GeneralMatrix*,MatrixType);
  98.  
  99. GeneralMatrix* AddedMatrix::Evaluate(MatrixType mt)
  100. { REPORT return GeneralAdd(bm1->Evaluate(), bm2->Evaluate(), mt); }
  101.  
  102. GeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mt)
  103. { REPORT return GeneralSub(bm1->Evaluate(), bm2->Evaluate(), mt); }
  104.  
  105. GeneralMatrix* MultipliedMatrix::Evaluate(MatrixType mt)
  106. {
  107.    REPORT
  108.    GeneralMatrix* gm2 = bm2->Evaluate();
  109.    MatrixType mtsym(MatrixType::Sym);
  110.    if (gm2->Type() == mtsym) gm2 = gm2->Evaluate(MatrixType::Rect);
  111.    return GeneralMult(bm1->Evaluate(), gm2, mt);
  112. }
  113.  
  114. GeneralMatrix* SolvedMatrix::Evaluate(MatrixType mt)
  115. { REPORT return GeneralSolv(bm1->Evaluate(), bm2->Evaluate(), mt); }
  116.  
  117. // routines for adding or subtracting matrices of identical storage structure
  118.  
  119. static void Add(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  120. {
  121.    REPORT
  122.    register real* s1=gm1->Store(); register real* s2=gm2->Store();
  123.    register real* s=gm->Store(); register int i=gm->Storage();
  124.    while (i--) *s++ = *s1++ + *s2++;
  125. }
  126.    
  127. static void Add(GeneralMatrix* gm, GeneralMatrix* gm2)
  128. {
  129.    REPORT
  130.    register real* s2=gm2->Store();
  131.    register real* s=gm->Store(); register int i=gm->Storage();
  132.    while (i--) *s++ += *s2++;
  133. }
  134.  
  135. static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  136. {
  137.    REPORT
  138.    register real* s1=gm1->Store(); register real* s2=gm2->Store();
  139.    register real* s=gm->Store(); register int i=gm->Storage();
  140.    while (i--) *s++ = *s1++ - *s2++;
  141. }
  142.  
  143. static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm2)
  144. {
  145.    REPORT
  146.    register real* s2=gm2->Store();
  147.    register real* s=gm->Store(); register int i=gm->Storage();
  148.    while (i--) *s++ -= *s2++;
  149. }
  150.  
  151. static void ReverseSubtract(GeneralMatrix* gm, GeneralMatrix* gm2)
  152. {
  153.    REPORT
  154.    register real* s2=gm2->Store();
  155.    register real* s=gm->Store(); register int i=gm->Storage();
  156.    while (i--) { *s = *s2++ - *s; s++; }
  157. }
  158.  
  159. static GeneralMatrix* GeneralAdd(GeneralMatrix* gm1, GeneralMatrix* gm2,
  160.    MatrixType mtx)
  161. {
  162.    int nr=gm1->Nrows(); int nc=gm1->Ncols();
  163.    if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
  164.       MatrixError("Incompatible dimensions");
  165.    if (!mtx) mtx = gm1->Type() + gm2->Type();
  166.    if (mtx == gm1->Type() && mtx == gm2->Type())
  167.                          // modify when band or sparse
  168.                          // matrices are included.
  169.    {
  170.       if (gm1->reuse()) { REPORT Add(gm1,gm2); gm2->tDelete(); return gm1; }
  171.       else if (gm2->reuse()) { REPORT Add(gm2,gm1); return gm2; }
  172.       else
  173.       {
  174.          REPORT
  175.          GeneralMatrix* gmx = gm1->Type().New(nr,nc); gmx->ReleaseAndDelete();
  176.          gmx->ReleaseAndDelete(); Add(gmx,gm1,gm2); return gmx;
  177.       }
  178.    }
  179.    else
  180.    {
  181.       if (gm1->Type()==mtx && gm1->reuse() )    // must have type test first
  182.       {
  183.      REPORT
  184.      MatrixRow mr1(gm1, StoreOnExit+LoadOnEntry+DirectPart);
  185.      MatrixRow mr2(gm2, LoadOnEntry);
  186.      while (nr--) { mr1.Add(mr2); mr1.Next(); mr2.Next(); }
  187.          gm2->tDelete(); return gm1;
  188.       }
  189.       else if (gm2->Type()==mtx && gm2->reuse() )
  190.       {
  191.          REPORT
  192.      MatrixRow mr1(gm1, LoadOnEntry);
  193.      MatrixRow mr2(gm2, StoreOnExit+LoadOnEntry+DirectPart);
  194.      while (nr--) { mr2.Add(mr1); mr1.Next(); mr2.Next(); }
  195.      if (gm1->Type()!=mtx) gm1->tDelete();
  196.          return gm2;
  197.       }
  198.       else
  199.       {
  200.          REPORT
  201.      GeneralMatrix* gmx = mtx.New(nr,nc);
  202.          MatrixRow mr1(gm1, LoadOnEntry);
  203.      MatrixRow mr2(gm2, LoadOnEntry);
  204.      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  205.      while (nr--)
  206.      { mrx.Add(mr1,mr2); mr1.Next(); mr2.Next(); mrx.Next(); }
  207.      if (gm1->Type()!=mtx) gm1->tDelete();
  208.      if (gm2->Type()!=mtx) gm2->tDelete();
  209.          gmx->ReleaseAndDelete(); return gmx;
  210.       }
  211.    }
  212. }
  213.  
  214.  
  215. static GeneralMatrix* GeneralSub(GeneralMatrix* gm1, GeneralMatrix* gm2,
  216.    MatrixType mtx)
  217. {
  218.    if (!mtx) mtx = gm1->Type() - gm2->Type();
  219.    int nr=gm1->Nrows(); int nc=gm1->Ncols();
  220.    if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
  221.       MatrixError("Incompatible dimensions");
  222.    if (mtx == gm1->Type() && mtx == gm2->Type())
  223.                          // modify when band or sparse
  224.                                              // matrices are included.
  225.    {
  226.       if (gm1->reuse())
  227.       { REPORT Subtract(gm1,gm2); gm2->tDelete(); return gm1; }
  228.       else if (gm2->reuse()) { REPORT ReverseSubtract(gm2,gm1); return gm2; }
  229.       else
  230.       {
  231.          REPORT
  232.      GeneralMatrix* gmx = gm1->Type().New(nr,nc); gmx->ReleaseAndDelete();
  233.          gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2); return gmx;
  234.       }
  235.    }
  236.    else
  237.    {
  238.       if (gm1->Type() == mtx && gm1->reuse() )
  239.       {
  240.          REPORT
  241.      MatrixRow mr1(gm1, LoadOnEntry+StoreOnExit+DirectPart);
  242.      MatrixRow mr2(gm2, LoadOnEntry);
  243.      while (nr--) { mr1.Sub(mr2); mr1.Next(); mr2.Next(); }
  244.          gm2->tDelete(); return gm1;
  245.       }
  246.       else if (gm2->Type() == mtx && gm2->reuse() )
  247.       {
  248.          REPORT
  249.          MatrixRow mr1(gm1, LoadOnEntry);
  250.      MatrixRow mr2(gm2, LoadOnEntry+StoreOnExit+DirectPart);
  251.      while (nr--) { mr2.RevSub(mr1); mr1.Next(); mr2.Next(); }
  252.      if (gm1->Type()!=mtx) gm1->tDe